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Abstract 

We discuss recent developments in the formulation of a Statistical 
Mechanics approach to non thermal systems, such as granular media. 
We review a few important numerical results on the assessment of 
Edwards' theory and, in particular, we apply these ideas to study a 
mean field model of a hard sphere binary mixture under gravity, which 
can be fully analytically investigated. As a consequence, we derive the 
rich phase diagram and predict the features of segregation patterns of 
the mixture. 



1 Introduction 

The issue we discuss here, which recently raised considerable interest [1, 2, 3, 
4, 5, 6], is the possibility to describe non thermal systems, such as granular 
media, by using concepts from Statistical Mechanics, as proposed by Edwards 
in 1989 [1]. 

Granular media [7] are dissipative systems made of macroscopic grains, 
such as powders or sands, where gravity effects are usually much stronger 
than thermal agitation, which is thus neglected (i.e., Tb a th = 0). As much 
as thermal systems, their macroscopic properties are characterized by a few 
control parameters and their macrostates correspond to a huge number of 
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Figure 1: A schematic picture of the system whole configuration space, flrot, 
and the subspace of Inherent States, flis, i-e., mechanically stable states. 

microstates. In granular media at rest, these microstates are mechanically 
stable, "frozen", arrangements of grains. They have been called "inherent 
states" [2], (see Fig.l), because of their connections to "inherent structures" 
[8] of supercooled liquids quenched at zero temperature, i.e., local minima of 
the potential energy in configuration space (which were introduced for their 
important role in glassy properties). 

In thermal systems the space of microstates is explored by the presence 
of a finite Tt, at h- In granular media at rest T^th = 0, but the effects of a 
finite thermal bath (Tb at h > 0) can be obtained by an external drive (like a 
"shake") of amplitude T > 0. 

The above considerations outline the possibility to develop a Statistical 
Mechanics approach to describe granular media. In the present paper, this 
approach is introduced in Sect. 2 and its check, via Monte Carlo simulations 
in the schematic framework of lattice models for granular media, is reviewed 
in Sect. 3. Sect. 4 is devoted to the analytical investigation of a mean field 
model of a binary hard sphere mixture under gravity for granular media, 
treated with the above approach. Finally, mean field predictions about the 
system phase diagram and its segregation phenomena are discussed in Sect. 
5. 
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Figure 2: Our lattice grain models are subject to a Monte Carlo dynamics 
made of "tap" sequences. A "tap" is a period of time, of length r (the tap 
duration), where they evolve at a finite value of the bath temperature, Tr 
(the tap amplitude); afterwards grains are frozen in one of their inherent 
states. 



2 Probability distribution 

Edwards [1] proposed a method to individuate the probability, P r , to find the 
system in its inherent state r, under the assumption that these mechanically 
stable states have the same a priori probability to occur. The knowledge of 
P r has the conceptual advantage to substitute time with ensemble averages, 
and thus allowing the description of the system properties by use of few basic 
theoretical concepts, as in thermodynamics. A possible approach to find P r 
is as follows [2] . P r is obtained as the maximum of the entropy, 

S = - J2 P r lnP r (1) 

r 

with the macroscopic constraint that the system energy, E = J2 r P r E r , is 
given. This assumption leads to the Gibbs result: 

P r oc e - Pcon f Er (2) 

where (3 con f is a Lagrange multiplier, called inverse configurational tempera- 
ture, enforcing the above constraint on the energy: 

ficonf q jp 

S conf = ]nQ IS {E) (3) 
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Here, Qjs{E) is the number of inherent states with energy E. Thus, sum- 
marizing, the system at rest has Tb at h = and T con f = ^ 0. 



2.4 




Figure 3: The time average of the energy, e, and (inset) its fluctuations, 
Ae 2 , recorded at stationarity during a tap dynamics, as a function of the tap 
amplitude, Tr, in the 3D lattice monodisperse hard sphere model. Different 
curves correspond to sequences of tap with different values of the duration 
of each single tap, r . 

These basic considerations, to be validated by experiments or simulations, 
settle a theoretical Statistical Mechanics framework to describe granular me- 
dia. Consider, for definiteness, a system of monodisperse hard spheres of 
mass m. In the system whole configuration space Qxot, we can write Ed- 
wards' generalized partition function as: 

where T~Chc is the hard core interaction between grains, mgH is the gravity 
contribution to the energy (H is particles height), and the factor IT (7) is 
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a projector on the inherent states space Qis (see Fig.l): if r G Qjs then 
Il r (7) = 1 else Il r (7) = 1 — 7. Usual Statistical Mechanics, where is 
identified with TJ, ai /j, is recovered for 7 = 0; Edwards' "Inherent States" 
Statistical Mechanics is obtained for 7 — > 1, and by definition (3~ x is called 

T C onf [ 1 j 2] . 




Figure 4: The time averages of energy fluctuations, Ae 2 , when plotted as 
function of the time average of energy, e, collapse on a single master function 
for all the different values of tap amplitude and duration, T r and r , shown 
in Fig. 3. Ensemble averages independently calculated over the distribution 
eq.(2), (Ae 2 ), go on the same curve when plotted as a function of (e) (the 
black empty circles). Time and Edwards ensemble averages coincide to a 
very good approximation. Data refer to a 3D monodisperse hard-sphere 
model under gravity, sketched in the inset. 
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3 Lattice hard sphere models for granular me- 
dia 



The above scenario has been checked [2] by Monte Carlo simulations, in a 
schematic model for granular media: a system of monodisperse hard-spheres 
(with diameter v^o) under gravity, where the centers of mass of grains are 
constrained to move on the sites of a cubic lattice (of spacing a = 1, see 
inset Fig. 4). Grains are subject to a dynamics made of a sequence of "taps" 
[10] (see Fig. 2): a single "tap" is a period of time, of length r (the tap 
duration), where they evolve at a finite value of the bath temperature, Tr 
(the tap amplitude); afterwards grains are frozen in one of their inherent 
states. 

We showed that under such a tap dynamics the systems reaches a sta- 
tionary state where time averages can be replaced, with good agreement, by 
ensemble averages over the above generalized Gibbs distribution eq.(2). 

Time averages over the tap dynamics at stationarity of the energy, e, and 
energy fluctuations, Ae 2 , as a functions of Tr, are plotted in Fig.3, where 
different curves correspond to tap sequences of different tap duration To- Ap- 
parently, e and Ae 2 are functions of both Tr and To, thus, for instance, Tr 
does not result to be a good parameter to uniquely individuate the stationary 
status of the system. Anyway, we can see that a single thermodynamic pa- 
rameter is necessary to characterize the macroscopic status. This is apparent 
from Fig. 4 (see also Fig. 5), where we plot Ae 2 as a function of e: all the 
data from different tap dynamics of Fig.3 collapse on a single master curve. 
We independently calculated the ensemble averages over distribution (2) of e 
and Ae 2 : interestingly these collapse on the same curve drawn by data from 
tap dynamics (see Fig. 4). 

Alternatively, an intensive "thermodynamic" parameter can be defined 
from the data collected along the tap dynamics at stationarity through the 
static fluctuation dissipation relation, which leads to individuate a "tem- 
perature", Tfd [2]. We numerically showed that Tfd coincides with Edwards' 
"configurational temperature", T con f, defined above. These results are briefly 
summarized in Fig. 5. There we show that averages of macroscopic quanti- 
ties calculated along different tap dynamics collapse on a single master curve 
when plotted as a function of the "thermodynamic" parameter Tfd- In Fig. 
5, we also show that ensemble averages, independently calculated by eq.(2), 
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collapse on the very same master curve when plotted as a function of T con f. 
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Figure 5: The density self-overlap function, Q, and (upper inset) the system 
density on the bottom layer, pt,, (previously recorded during tap dynamics 
with different durations tq and amplitude Tr) are plotted as function of 
Tfd (in units mgao). These time averages are compared with the ensemble 
averages over the distribution eq.(2) (the black empty circles), plotted as 
function of T con f (in units mga ). There is a very good agreement between 
the independently calculated time averages over the tap dynamics and the 
statistical mechanics ensemble averages a la Edwards. Data refer to a 3D 
monodisperse hard-sphere model under gravity. 

The tap "dynamics" is interesting in itself [10, 3, 15]. For instance, density 
relaxation is well fitted by stretched exponential and by logarithmic laws 
respectively at high and low Tr [2, 10], a fact also experimentally found [11]. 
In particular, one observes that by lowering Tr the systems gets "jammed" , in 
the sense that its characteristic relaxation time, r, enormously grows. This is 
shown in Fig.6. The presence of strong "aging" effects in models for granular 
media was indeed discovered a few years ago [3, 15]. The growth of r with Tr 
is well approximated by an Arrhenius or Vogel-Tamman-Fulcher law [2, 10], 
resembling the slowing down of glass formers close to the glass transition, a 
result recently experimentally confirmed [12]. 

Eq.(4) gives a general framework where the properties of thermal and 
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Figure 6: Inset The density self-overlap, Q, as a function of the number of 
taps, t n , for several Tp. Main Frame The characteristic relaxation time, r, 
defined by a stretched exponential fit of the long time dynamics of Q(t n ). At 
low Tr, r appear to diverge a la Arrhenius or Vogel-Tamman-Fulcher (the 
superimposed fit). 

non-thermal systems, such as their phase diagram, can be discussed. For 
instance, a relation, as the one sketched in Fig. 7, is now found to exist 
between glassy properties of hard spheres under gravity in usual statistical 
mechanics (7 = limit) and the "jamming" properties of granular media (7 = 
1 limit); the control parameters being the number of grains and, respectively 
the temperature, Tb at h, and the configurational temperature, T con f. This is 
consistent with the theoretical speculation by Liu and Nagel [9] about the 
existence of a generalized phase diagram for the "jammed" region of granular 
media and "glassy" region of glass formers. 

Also in glasses, following for example the inherent structure approach 
[8], one can define a configurational entropy associated to the number of 
inherent structures corresponding to a fixed energy, E, and consequently the 
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configurational temperature. When the system is frozen at zero temperature 
in one of its inherent states it does not evolve anymore. A way to visit the 
inherent structures is by using a tap dynamics (i.e., a procedure similar to 
that used above for granular materials). Thus, more generally, the framework 
of eq. (4) opens the way to a Statistical Mechanics of "frozen system" . 

Finally, it is important to recall that we have shown [2] that more than 
one Lagrange multiplier may be necessary, in general, to assign the macro- 
scopic system status. For instance, in the case of a binary mixture of hard 
spheres on a lattice, at least two configurational temperatures must be intro- 
duced. In this case, an extension of Edwards' approach is required, where 
P r is given by the maximum of the entropy with separate constraints on the 
gravitational energies of the two species of grains: E\ = m\g P T H\ r and 
E 2 = rn 2 g^2 r P r H 2r (Hi and H 2 are the heights of species 1 and 2). This 
gives two Lagrange multipliers: 

. _ dlnn IS (E 1 ,E 2 ) dlnQ IS (E 1 ,E 2 ) 

Pi — KT=, Pi 



dE 1 ^ dE 2 

where fi/s(Ei, E 2 ) is the number of inherent states with Ei,E 2 . Edwards' 
generalized partition function must now be written as: 

% = e -KHC-l3im 1 gH 1 -l32m2gH 2 . ^(7) 

To illustrate the theoretical developments allowed by Edwards formula- 
tion, we now discuss a mean field model of the above binary mixture. We can 
thus predict the structure of its phase diagram and describe an ubiquitous 
and intriguing phenomenon, called "size segregation", observed in vibrated 
granular mixtures where in presence of shaking the system is not randomized, 
but its components tend to separate. 



4 A Mean Field model 

We consider a lattice of H layers, occupied by a hard spheres system made 
of two species, 1 (small) and 2 (large), with grain diameters ao and V2oq 
(the lattice spacing is ao = 1). The hard core potential Tine is such that 
two nearest neighbor particles cannot overlap. This implies that only pairs 
of small particles can be nearest neighbors on a layer. 
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Figure 7: A schematic picture of a unifying phase diagram for the "glassy 
region" of ordinary hard spheres under gravity and "jammed region" of gran- 
ular media which emerges from the Statistical Mechanics approach to pow- 
ders. 



Each layer is dealt with the standard mean field technique, i.e., it is 
subdivided in two sublattices A and B, where each site in sublattice A (resp. 
B) interacts with all sites in B (resp. A) and only with them (see Fig.8). 
The Hard-Core lattice Hamiltonian for spheres on the same layer at height 
z (with z G {1, H}) is thus: 

n B o{z) = ± y £f{ n t',nf>) (5) 

ij 

Here we introduce the notation nf z (resp. nf z ) to individuate the status of 
site i in sublattice A (resp. B) of the layer at hight z: nf' B G {0, 1, 2} is an 
occupancy variable which is if no grains are on site i, 1 or 2 whether a small 
or a large sphere is there. In eq.(5) the sum is on all sites of sublattices A z 
and B z , i,j G {1, ...iV}. The "shape factor", f(rii,rij), indicates whether two 
grains can be nearest neighbor: /(n^n,) = (1 - S nit0 )(l - S njt0 ) - S ni!l 8 njt i. 
The Hard-Core limit is found for J — ► oo. 

To set the average density a chemical potential term is also present: 

ft-ni*) ^ 0Mnf,i + ^2^,2) (6) 

i,a 
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Here the sum is on all sites of sublattices A z and B z (i.e, % G {1, ...N} and 
a G {A Z ,B Z }). Of course, also the gravity contributions are considered: 



rriigz 



(7) 



i. a 



where z is the height of the layer where belong to, and m^, k G {1,2}, is 
the mass of grains of type 1 or 2. 




Figure 8: A schematic picture of the lattice considered here. We have H 
horizontal layers, each divided in two reciprocally fully connected sublattices 
A and B. Grains are located on the lattice sites and, in order to be "stable" , 
must rest on their homologous underneath grains. 



4.1 Edwards projector term 7iEdw 

The above Hamiltonian is the usual one describing a set of hard spheres un- 
der gravity at a mean field level. The novelty we consider here is Edwards' 
idea that in granular media averages must be done by considering only me- 
chanically stable states (i.e., "Inherent States", IS). In the present case of 
a lattice model the "stability" of a grain is given by the presence of some 
underlying grains which support it. The projector on IS has thus a simple 
form: it can be written as a sum of single particle terms to be added in the 
Hamiltonian. In fact, to decide whether a grain on layer z is "mechanically 
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stable" we must just check the status of its corresponding sites in the two 
layers below it. In the case of not fully connected layers also next neighbors 
sites at level z — 1 should be considered, but this is not the case here. In 
particular, since for a matter of simplicity we just discuss a system made of 
2 layers (i.e., H = 2), only the status of site i at z — 1 must be considered to 
establish whether at z is stable or not. The case with more than 2 layers 
gives very similar results to those presented here, but at the price of a more 
intricate analytical treatment. Since our model is already very schematic, 
we disregards these details and restrict here to H = 2. Thus, summarizing, 
in the present case the "stability" of grains on level z is given by the follow- 
ing term, which couples adjacent layers between them, to be added to the 
Hamiltonian [13]: 

H Edw (z) = Kj^i 1 - *n?«,oK».-i i0 (8) 

i,a 

Usual Statistical Mechanics is recovered for K = 0. Edwards' "Inherent 
States" Statistical Mechanics is obtained for K — > oo (i.e., K = ln(l— 7)); the 

projector on Inherent States being: IT({n" z }, K) = exp | — Ylf=i ^Edw( z )^- 
The system partition function can thus be written as: 

Z = £expj-^[Wtf C (*)+/3i?4 1 >(*)+^ 

{nf z } I 2=1 



5 Phase Diagram and Segregation Phenom- 
ena 

After standard manipulation of eq.(9), the saddle point free energy is ob- 
tained: 



z=l Lj.fc=l>2 



X A \ Z ) X B \Z) ~ X A K Z ) X B K Z ) 



E 

,fc=i,: 

-\n(Z[x A )Z[x B ]) (10) 
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where, after the K — > oo limit, we have defined 

Z[x A ] = 1 + e~ J[x( *W ] e Ml + e ~ J[ *A )(1)+!B 2 ) ( 1 )] e^ 2 + Q-J^Tw+^Ti 2 )] e 2i*x-Pimigho 

_|_ e -J[^ ) (l)+^ ) (l)+x^(2)] e tM 1 +ti 2 -Pim 1 gh + e -J[x^ ) (l)+x^ ) (2)+x ( 2 ) (2)] ^x+^-fhm^gho 
+ e -J[4 1) (l)+4 2) (l)+4 1) (2)+4 2) (2)] tflto-frmagho ( U ) 

Here we have introduced the density field a?i^(z), which is the number density 
fraction of species k G {1, 2} on sublattice a G {A, B} at level z G {1, if}. 
The constant ho is the vertical layer spacing. The 4H mean field self con- 
sistent coupled equations, to be solved in the limit J — > oo, are the saddle 
points: 

df[x A ,x B ] 
dx^ (z) 

Without entering further details, eq.s(12) admit three kinds of solutions, 
corresponding to three different possible phases for the system. Phase I 
(which is twofold) is characterized by the property that there are no large 

(21 (2) 

grains sitting on the bottom layer, i.e., x A '(1) = 2^(1) = 0. Phase II 
(twofold) admits the presence on both layers of both kind of particles. Finally, 
Phase III has the property that small grains have expelled large grains, i.e., 
x^a\z) = Va, z. The resulting phase diagram of the system is shown in 
Fig. 9. In general the crossing from one to the other phase is marked by 
a first order phase transition. Finally, it is possible to show [14] that the 
"Bethe lattice" version of the model used in the present "fully connected" 
mean field approach appears to be able to locate also a glassy phase in the 
above phase diagram, as originally predicted in [10, 15, 3]. 

It is worth to exhibit the configurational entropy of the system: 

Sconf = [fiidp! + fodfo + A*i<9 Ml + ^ 2 <9 M2 - 1] / (13) 

which is shown in Fig. 10. Finally, we discuss "size segregation", and specif- 
ically, the phenomenon known as the "Brazil Nut Effect" (BNE), i.e., the 
tendency of large grains to segregate on top of the system. The vertical 
segregation parameter is the difference of average heights of species 1 and 2: 
Ah = hi — It is possible to show that in Phase I large grains are always 
on average on top, i.e., Ah < 0, due to the geometric organization of the 



(12) 
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pack. In Phase II, BNE is observed when miPi > "^12^2-, else the "Reverse 
Brazil Nut Effect" (RBNE) is found (where large grains are on the bottom). 
Notice that in Phase II a route from the BNE to the RBNE region does not 
necessarily cut any transition line. Instead, crossing from the BNE region 
of Phase II to BNE region of Phase I (as in panel a) of Fig. 9) corresponds 
to cross a 1st order phase transition line (with a jump in Ah, as shown in 
Fig. 11). Phase III is associated to an other form of segregation where, as 
already stated, small grains fully expel large grains. This phase separation 
phenomenon can give rise to "vertical" as well as "horizontal segregation", 
i.e., pattern formation. 

Some numerical simulations [2, 16, 18] and experimental observations 
[17] about BNE appear to be explained by segregation mechanisms as those 
found in Phase I, II and III of the present model. In particular, segregation 
under "condensation", numerically discovered in [18], may be originated by 
phenomena as those found in Phase III. 

6 Conclusions 

As a theoretical picture is taking shape, there are still many open problems 
with the Statistical Mechanics of granular media; of course, the assumption 
that the stationary distribution is given by the maximum of an entropy with 
a few macroscopic constraints is to be checked case by case; and there is 
no general a priori criteria to select such necessary macroscopic constraints, 
i.e., there is no general criteria to establish how many "thermodynamic" 
parameters (T con /s) are needed to characterize the status of a given powder 
[2, 19]. 

Anyway, simple models of granular media have confirmed and extended 
[2, 3, 4, 5, 6] Edwards' 1989 hypothesis [1] that granular media in their "in- 
herent states" can be described by a generalized Statistical Mechanics. Here 
we have reviewed recent progresses made in this research area and, in par- 
ticular, we have analytically investigated a simple mean field model where 
to apply Edwards' Statistical Mechanics approach. A comprehensive picture 
emerges of "thermodynamical" and dynamical properties of these non ther- 
mal systems, ranging from segregation patterns to their jamming transition 
[1, 2, 20]. The presence of phase transitions in the inherent states space, here 
predicted by a mean field theory, is to be confirmed by experimental work. 
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Figure 9: Phase diagrams of the mean field hard spheres mixture model in 
the plane (e - ™ 1 ^ 1 , e~ m2 ^ 2 ) described in the text. In panel a), b), c), d) the 
parameters (e^e^ 2 ) are respectively: (.1, .1); (1, .1); (.1, 1); (1, 1). In Phase I 
only BNE is present for any value of mi/3i and m^fti- in Phase II crossover 
from BNE to RBNE occurs when mi/3i/m 2 /52 = 1 (the diagonal in the graph). 
In Phase III large grains are expelled from the system. 
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Figure 10: The configurational entropy, s con f, of the mean field model as a 
function of e" mi/31 for (e^,e^ 2 ) = (1, 1) and e"™ 2 * = .5. 
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Figure 11: The vertical segregation parameter, Ah = h\—h2 (i.e., the average 
height difference of species 1 and 2) of the mean field model, as a function 
of e -mi()i for ( e w ; e^ 2 ) = (.l,.l). Left panel: e"™ 2 * = .75; right panel: 
e -m 2 /? 2 = _ 5g _ 
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